Background

We looked for combination therapies of neglected and misused antibiotics. Rifampin and polymyxin B combination stood out as a promising approach against scpectrum of clinical isolates. This repository contains all the data and code to reproduce the analysis. For details and full captions of figures, consult the publication.

Combinations

We studied the potential of the rifampicin-polymyxin B combination against intra- and extracellular forms of bacteria: three P. aeruginosa strains, two clinical isolates of A. baumannii, E. coli, and K. pneumoniae.

dat = loadDdiData('input/dat/rds/ddi_all_strains.rds')

# Response surface analysis
# rs = getRS(dat, 'rs')
# rs = readRDS('input/dat/tmp/rs.rds')
# 
# open3d()
# lapply(names(rs), sav3D)

#getDdiCounts(readRDS('input/dat/tmp/rs.rds')) %>% 
#  group_by(Strain) %>% 
#  ggplot(aes(Cond, Syn)) +
#  geom_col(fill = 'grey70', col = 1) +
#  facet_wrap(~Name, nrow = 3, dir = 'v') +
#  labs(y = '# of synergistic concentrations', x = '') +
#  theme(strip.text.x = element_text(size = 12)) 
#ggsave(filename='output/fig/SynergyCountAcid.svg', bg='white', h=7.5, w=7.5)

Dose Response And pH

theme_set(theme_pdf())

sub = subset(dat, strain == 'ATCC27853' & (d1 == 0 | d2 == 0))

p1 = function() pltDR(effect ~ c1, subset(sub, c2==0), xlab = 'RIF, mg/L', zero=30)
p2 = function() pltDR(effect ~ c2, subset(sub, c1==0), col='orange', xlab = 'PMB, mg/L')

# pH
foo = fread("input/dat/raw/growth_and_pH.tsv")

foo = foo[time_h >= 0 & overnight_medium == 'pH7.4', .(
  od = gmu(od620),
  od.ci = gci(od620),
  ph = mean(pH),
  ph.ci = ci(pH)
  ), .(growth_medium, time_h)]

limit_of_detection = 0.03
foo[, od := replace(od, od < limit_of_detection, limit_of_detection)]

p0 = foo %>% 
  mutate(
    growth_medium = relevel(factor(growth_medium), 'pH7.4'),
    time_h = ifelse(time_h > 20, 20, time_h),
    NULL
    ) %>% 
  ggplot(aes(time_h, od, shape = growth_medium)) +
  scale_x_continuous(breaks = c(0, 5, 10, 20), labels = c(seq(0, 10, 5), 24), limits = c(0, 20)) +
  scale_shape_manual(name = '', values = c(19, 22, 15), labels = c('pH 7.4', 'pH 5.5',
      'pH 5.5 buffered')) +
  theme(legend.position = 'top')

p3 = p0 + 
  scale_y_continuous(trans='log2', breaks = c(0.01, 0.03, 0.1, 0.3, 1, 3,
      10)) + 
  geom_point(size=.point_size, col='grey20') +
  geom_pointrange(aes(ymin = od*od.ci, ymax = od/od.ci)) +
  geom_line(col='grey20') +
  coord_cartesian(ylim = c(0.01, 10)) +
  labs(x = 'Time, h', y = expression('OD'[620]))

p4 = p0 %+% aes(y = ph) + 
  scale_y_continuous(position = 'right') +
  geom_point(size=.point_size, col=cbPalette[3]) +
  geom_line(col=cbPalette[3]) +
  geom_pointrange(aes(ymin = ph-ph.ci, ymax = ph + ph.ci), col=cbPalette[3]) +
  theme(legend.position = 'none',
    axis.title.y.right=element_text(color = cbPalette[3]),
    axis.text.y.right=element_text(color = cbPalette[3]),
    axis.ticks.y.right=element_line(color = cbPalette[3])
  ) +
  labs(x = '', y = 'pH')

foo = fread('input/dat/raw/charcoal.csv')
labs = c('RIF, mg/L', 'PMB, mg/L', 
  rep(expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')'), 2)
)

tmpFun = function (i, df) {
  par(mar=c(3.1, 6.1, 0, 0.7), mgp = c(1.7, 0.5, 0), tcl=-0.2)
  if(i == 1) bar = subset(df, ab == 'rif')
  else bar = subset(df, ab == 'pmb' | (ab == 'rif' & conc == 0))

  # fit char vs no-char but share Emin and EC50 (argument pmodels)
  fit = drm(response~conc, data=bar, CURVE, 
    pmodels=data.frame(CURVE,1,1),  
    fct=LL2.4( fixed=c(1, NA, NA, NA))) 
  
  plot(fit, type="average", cex.axis=1.25, cex.lab=1.25, lwd=2, 
    col='grey30', pch = c(1, 19), lty = 1, asp=1,
    xlab = labs[i], ylab = labs[i+2], broken=T,
    ylim=c(-4,4), legendText=c("control", "charcoal")
  )
  plot(fit, type="confidence", add=T, legend=F, 
    col=c('grey70', 'grey10')
  )
}

p5 = function() tmpFun(1, foo)
p6 = function() tmpFun(2, foo)

aligned = align_plots(p3, p4, align = 'hv', axis = 'tblr')

plot_grid(
  ggdraw(p1), ggdraw(p2), ggdraw(aligned[[1]]) + draw_plot(aligned[[2]]),
  plot_grid(ggdraw(p5), ggdraw(p6), nrow=2),
  labels = 'AUTO', align = 'v'
)
# ggsave('output/fig/SFig_mth_dr_pH_charcoal.svg', width=10, height=9)

PD parameters

# schematic for PD pars
tmp = data.table(
  #Conc      = c(0, 0.1, 0.3, 1, 3, 10, 30, 100),
  Conc      = c(0, 0.01, 0.03, 0.1, 0.3, 1, 3, 10),
  Response  = c(3, 3, 2.8, 1.5, -0.5, -2.5, -2.8, -3.2)
)

tmp[, Response := Response + 1]

# svg('output/fig/SFig_PDpar.svg', width = 3.5, height = 3.5)
#oldpar = par()
par(tcl = -0.3, mar=c(4.1, 4.1, 1.5, 0.5), mgp=c(1.5, 0.75, 0))

plot(
  (mod = drm(Response ~ Conc, data = tmp, 
      fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
      lowerl = c(-2.5, NA, NA),
      robust = 'median'
  )),
    broken = T, pch = 19, 
    col = 'grey40',
    ylim = c(-2.5, 4.7), xlim = c(0, 50),
    xlab = c('Conc, mg/L'),
    ylab = expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')'),
    asp = 1
)

x1 = ED(mod, 0, type = 'absolute', logBase = exp(1))[1]
## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:0 0.4249298  0.0092374
x2 = ED(mod, 50, logBase = exp(1))[1]
## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:50 0.233915   0.005085
y2 = predict(mod, data.frame(Conc = x2))

segments(x1, 0, x1, -100, lty = 2)
segments(0.0001, 0, x1, 0, lty = 2)
segments(x2, y2, x2, -100, lty = 2)
segments(0.0001, y2, x2, y2, lty = 2)

abline(h = coef(mod)[2], lty = 2)
abline(h = coef(mod)[1], lty = 2)

text(x1+0.3, y = 0 + 0.1, expression(C[s]))
text(x2+0.3, y = y2 + 0.2, expression(EC[50]))
text(x = 0.002, coef(mod)[2] + 0.3, expression(E[min]))
text(x = 0.002, coef(mod)[1] + 0.4, expression(E[max]))
# dev.off()
# add Tn screen's PA14 (Liberati), and pmbR strains PA947 and PA1292 
foo = rbind(
  # pa14
  fread('input/dat/raw/pa14_mut_pmb_rif_cfu.csv')[mut == 0] %>% 
    .[, `:=` (genus   = 'Pseudomonas',
              strain  = 'PA14 Liberati',
              mut     = NULL)],
  # pmbR
  fread('input/dat/raw/pmb_resistant_isolates.csv') %>% 
    .[, genus := 'Pseudomonas']
  ) %>% 
  .[, `:=` (
    c1 = d1*mic1,
    c2 = d2*mic2)]

# transform
dat_fit = rbind(foo, dat, fill = T)[c1 == 0 | c2 == 0,
  ][, dose := ifelse(d1 == 0, c2, c1)
  ][, ab   := ifelse(d1 == 0, 2, 1)
  ][, .(data = list(.SD)), .(strain, cond)]

# fit four parameter logistic regression
dat_fit[, fit := map(data, ~drm(effect ~ dose, ab, data = .x, 
  fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
  pmodels = data.frame(ab, 1, ab),
  lowerl = c(-4.5, NA, NA)))]

# get parameters
dat_fit[, pars  := map(fit, ~getPDPar(.x))]
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:2:0    1.0936    0.7459    1.6033
## e:1:0 3017.5829 1432.1220 6358.2618
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:1:0 199.80084 143.78167 277.64579
## e:2:0   1.09877   0.85522   1.41167
## 
## Estimated effective doses
## 
##         Estimate      Lower      Upper
## e:1:0  4195.4112   844.9468 20831.4600
## e:2:0     1.4129     1.1031     1.8096
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  136.099  94.119 196.802
## e:2:0  108.833  77.678 152.485
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  116.627  78.836 172.534
## e:2:0  100.888  71.106 143.146
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  2.19824 1.73299 2.78839
## e:2:0  1.05512 0.74755 1.48923
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0  0.54794 0.41895 0.71664
## e:1:0  2.54461 2.01152 3.21898
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  1.23558  0.90005  1.69619
## e:1:0 12.27513  9.50191 15.85774
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  0.72807  0.50460  1.05050
## e:1:0  8.19634  6.45160 10.41292
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0   8.8108  5.2898 14.6755
## e:1:0  29.9697 20.8806 43.0152
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0   2.9223  1.9854  4.3015
## e:1:0  25.5587 19.0791 34.2388
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  1.08515  0.83477  1.41063
## e:1:0 17.68410 13.82126 22.62653
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0   3.0869  1.4710  6.4777
## e:1:0   9.1434  7.0217 11.9061
## 
## Estimated effective doses
## 
##       Estimate  Lower  Upper
## e:2:0   2.7579 1.6050 4.7391
## e:1:0   6.2385 4.5268 8.5975
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  0.71158 0.43400 1.16668
## e:2:0  4.78692 3.03655 7.54625
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0  0.61538  0.37753  1.00308
## e:2:0 11.63563  7.65409 17.68830
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 10.31467  8.03876 13.23493
## e:2:0  0.52524  0.40478  0.68156
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 21.90092 17.02990 28.16518
## e:2:0  0.58875  0.45974  0.75395
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0   34.814  26.015  46.590
## e:2:0  212.935 151.138 300.000
## 
## Estimated effective doses
## 
##       Estimate  Lower  Upper
## e:1:0   55.495 42.649 72.210
## e:2:0   24.544 18.867 31.930
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  22.7675 16.2560 31.8872
## e:2:0   2.0841  1.4483  2.9992
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 18.17315 12.63799 26.13261
## e:2:0  0.36616  0.23469  0.57127
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:1:0  6.273396  3.311455 11.884652
## e:2:0  0.071390  0.039096  0.130359
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  3.89931 2.70443 5.62209
## e:2:0  0.50220 0.38022 0.66331
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  5.68227 3.61600 8.92927
## e:2:0  1.27999 0.91009 1.80022
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  15.6419 11.5468 21.1894
## e:2:0   1.1580  0.8129  1.6495
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 34.12972 26.34817 44.20944
## e:2:0  0.56378  0.45798  0.69403
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  63.2718 45.7219 87.5581
## e:2:0   6.4509  4.8977  8.4967
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:1:0  78.80048  58.25397 106.59387
## e:2:0   1.15853   0.89306   1.50291
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  19.0752 14.1575 25.7011
## e:2:0   3.9380  2.8444  5.4521
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 38.15049 27.24067 53.42967
## e:2:0  0.92200  0.66327  1.28166
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  0.60558  0.38290  0.95774
## e:1:0 27.16329 14.77923 49.92439
dat_fit = dat_fit[
  dat_fit[, data[[1]][1, .(genus, mic1, mic2)], .(strain, cond)], 
  on = c('strain', 'cond')]

# reformat
dat_fit = unnest(dat_fit, pars) %>% as.data.table
dat_fit = getPdReadyForPlotting(dat_fit, dat)

p1 = dat_fit[!grepl('C', par)] %>% 
  .[lwr < -10, c('lwr', 'upr', 'value') := NA] %>% 
  ggplot(aes(value, strain, color = genus)) +
    geom_vline(xintercept = 0, col = 'white') +
    geom_pointrange(aes(xmin = lwr, xmax = upr), size = 1) +
    geom_text(aes(label = round(signif(value, 2), 1)), col = 1, size = 2.5) +
    facet_grid(cond ~ par, labeller = label_parsed, 
      scales = 'free', space = 'free_y') +
    geom_blank(aes(x = min)) +
    geom_blank(aes(x = max)) +
    scale_color_manual(values = cbPalette[4:1], guide = guide_legend(
        label.theme = element_text(angle = 0, face = "italic"), reverse = T)) +
    scale_fill_manual(values = cbPalette[4:1]) +
    labs(y = '', color = '') +
    theme_html() +
    theme(legend.position = 'top')

p2 = p1 %+%
  mutate(dat_fit[grepl('EC', par)], 
    value = ifelse(value > 512, 512, value), 
    mic = ifelse(mic > 512, 512, mic)) +
  geom_point(aes(mic, col = genus), size = 4.5, pch = 0) +
  geom_text(aes(x = mic, label = round(mic, 1)), size = 2.5) +
  scale_x_continuous(trans = 'log10', limits = c(0.03, 512), 
    labels = function(x) ifelse(x == 0, "0", x)) +
  facet_grid(cond ~ par, labeller = label_parsed, 
    scales = 'free', space = 'free_y') +
  xlab('Conc, mg/L') +
  theme(legend.position = 'none')

plot_grid(p1 + xlab(expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')')), 
  p2, rel_widths = c(3, 2), align = 'hv', axis = 'tb')
# ggsave('output/fig/SFig_mth_pd.svg', width=10, height=6)

Time-kill

dtk = fread('input/dat/raw/time-kill.csv') %>% 
  remNonCharcoal() %>% 
  addNewVariablesToTimeKill()

# pltTimeKill(dtk)
#ggsave('output/fig/SFig_TimeKill.svg', width=15, height=5)

# # Print out the concentrations used
# unique(dtk[time_h > 0], by = c('broth', 'group')) %>% 
#   .[, .(broth, group, d_1, d_2)] %>% 
#   .[order(group)]

Reverse Genetics Screen

To account for the synergy in molecular terms—beyond a nonspecific increase in membrane permeability by polymyxin B—we turned to chemical genetics (Brochado and Typas, 2013). Working with ordered PA14 transposon library (Liberati et al., 2006), we derived a growth measure for monotherapies and combinations using colony opacity (Kritikos et al., 2017). To account for plate-to-plate variation, the opacity was multiplicatively corrected. This results in zero-centering of the Bliss scores, which were derived next. The significance of difference from zero Bliss score, for any mutant, was estimated by a T-test (5 biological replicates) and corrected for multiple testing (Benjamini-Hochberg).

tmp = chem.gen.res %>% 
  mutate(t.test.q.value = ifelse(t.test.q.value < 0.001, 0.001,
      t.test.q.value)) %>% 
  mutate(showname = ifelse(
      (t.test.q.value < 0.05 & (gene.name != '' | gene.name.to.show %in% c('66480', '26590', '43270', '02150'))) | t.test.q.value < 0.01 ,
      T, F)) 

# tmp[t.test.q.value < 0.05 | t.test.q.value.gene < 0.05, .N, media]
# #      media  N
# # 1:      LB 65
# # 2: LBpH5.5 39

# # only 3 are shared between the media
# foo = tmp[media == 'LB' & (t.test.q.value < 0.05 | t.test.q.value.gene < 0.05), Tn.mutant.id] 
# bar = tmp[media == 'LBpH5.5' & (t.test.q.value < 0.05 | t.test.q.value.gene < 0.05), Tn.mutant.id] 
# tmp[Tn.mutant.id %in% intersect(foo, bar)]


getLocus = function (dat, condition) {
  unique(dat[eval(parse(text=condition))], by = c('locus'))[, locus] 
}

by_gene       = getLocus(tmp, 't.test.q.value.gene < 0.05')
by_mut        = getLocus(tmp, 't.test.q.value < 0.05')
by_both       = getLocus(tmp, 't.test.q.value < 0.05 | t.test.q.value.gene < 0.05')
by_gene_only  = setdiff(by_gene, by_mut)

# # 70 hits
# length(by_both)
# # 53 from mutants
# intersect(by_mut, by_both) %>% length()
# # 17 additional from genes
# setdiff(by_gene, by_mut) %>% length() 
# # 31 shared between mutants and genes 
# intersect(by_gene, by_mut) %>% length()

Of the 70 hits, 53 are from mutant level of analysis and 17 come from gene level of analysis. 31 are shared between mutants and genes i.e. hits regardless if mutant or gene level of analysis is used.

(p0 = filter(tmp, media == 'LB') %>%  
  {
  ggplot(., aes(epsilon.median.mutant.condition, -log10(t.test.q.value), col = t.test.q.value < 0.05)) +
  geom_point(size = 2.5) +
  scale_color_manual(values = c('grey70', 'grey20')) +
  ggrepel::geom_text_repel(data= . %>% .[(showname)], 
    aes(label=gene.name.to.show), size = 6,
    max.overlaps = Inf, box.padding=0.75, seed=3) +
  lims(x = c(-0.8, 0.6), y = c(0, 3)) +
  labs(x = 'Shift in Bliss score', y = expression('log'[10]*'(p-value)')) +
  theme(aspect.ratio = 1, legend.position = 'none')

})
# ggsave('output/fig/ChemGen_LB.svg', width=5, height=5)
p0 %+% filter(tmp, media == 'LBpH5.5')
# ggsave('output/fig/ChemGen_LBpH5.5.svg', width=5, height=5)
LB LB pH 5.5
getGenes = function (dat) {
  out = list()
  for(i in unique(dat$media))
    out[[i]] = unique(dat[t.test.q.value < 0.05 & media == i], by='locus') %>%
      .[, (gene.name.to.show)]
  return(out)
}

foo = getGenes(chem.gen.res)
#vtable = gplots::venn(foo)

GO terms

With the following, we bring some biological knowledge into the analysis. This will get us at the level of processes/compartments as opposed to individual genes. We will focus on LB beacuse we have more data from there which is also more reliable (75% of the unique hits come from LB; there is less variance). In addition, our results suggest, the effect of pH on synergy is weak in PA14.

We did Gene Set Enrichment Analysis (GSEA) using GO terms from pseudomonas.com website and Kologorov-Smirnov testing for statistical significance estimation. Although most common approach, it has been critizised for example here.

Cell component

library('topGO')
tab = getTopGO(chem.gen.res)

tmpFun = function (dat, name='CC') {
  dat[[name]][1:10, c('GO.ID', 'Term', 'raw.p.value')] %>% 
    mutate(Pvalue = round(as.numeric(raw.p.value), 4))  %>% 
    dplyr::select(-raw.p.value) %>% 
    mutate(`GO.ID` = sub('GO:', '', `GO.ID`)) 
}

tmpFun(tab) %>%
  knitr::kable()
GO.ID Term Pvalue
0008076 voltage-gated potassium channel complex 0.0096
0016021 integral component of membrane 0.0132
0055052 ATP-binding cassette (ABC) transporter complex, substrate-binding subunit-containing 0.0193
0005694 chromosome 0.0362
0033573 high-affinity iron permease complex 0.0530
0009289 pilus 0.1160
0005839 proteasome core complex 0.1221
0005960 glycine cleavage complex 0.1396
0030257 type III protein secretion system complex 0.1701
0005615 extracellular space 0.1717
  # knitr::kable(format='latex', booktabs=T) %>% 
  # kableExtra::save_kable('output/fig/table_cell_component.png', font='Verdana')

Biological process

# biological processes
tmpFun(tab, 'BP') %>% 
  knitr::kable()
GO.ID Term Pvalue
0009236 cobalamin biosynthetic process 0.0086
0009116 nucleoside metabolic process 0.0090
1900727 osmoregulated periplasmic glucan biosynthetic process 0.0166
0019700 organic phosphonate catabolic process 0.0173
0070475 rRNA base methylation 0.0179
0044249 cellular biosynthetic process 0.0260
0019354 siroheme biosynthetic process 0.0298
0006437 tyrosyl-tRNA aminoacylation 0.0346
0044262 cellular carbohydrate metabolic process 0.0357
0019748 secondary metabolic process 0.0370
  # knitr::kable(format='latex', booktabs=T) %>% 
  # kableExtra::save_kable('output/fig/table_biological_process.png', font='Verdana')

Protein-protein interaction network

Protein-protein interaction (PPI) analysis using STRING database. There was no data on PA14, so we will use PAO1 data to build and analyse the network onto which we then map PA14 orthologs.

library(igraph)
library(stringr)

hits = addAnnotation(chem.gen.res[t.test.q.value < 0.05, .(colony)])
hit_list = unique(hits[PAO1.ortholog != '', PAO1.ortholog])

# bring in the interactome
string_11_pa = 'input/dat/raw/pao1_ppi_stringDB_287.protein.links.v11.0.txt' %>% 
  fread() %>% 
  .[, `:=` (
    combined_score = combined_score/1000,
    #strip / trim species ID
    protein1 = gsub("287.DR97_", "PA", protein1, fixed=TRUE),
    protein2 = gsub("287.DR97_", "PA", protein2, fixed=TRUE)
  )]

# create an igraph object
pa_ppi = graph.data.frame(string_11_pa, directed = FALSE)
# summary(pa_ppi)
V(pa_ppi)$degree = degree(pa_ppi)

pa_ppi = simplify(pa_ppi, remove.multiple=T, remove.loops=T, edge.attr.comb='mean')

# extract network for all the hits in the provided list
hit_ppi = induced.subgraph(pa_ppi, which(V(pa_ppi)$name %in% hit_list))

# store degree on the network (main hit network)
V(hit_ppi)$degree = degree(hit_ppi)

# Replace the PAO1 genes with corresponding PA14 ones
lookup = unique(hits[, .(PAO1.ortholog, gene=gene.name.to.show)], by='PAO1.ortholog')
V(hit_ppi)$name = left_join(data.table(PAO1.ortholog=V(hit_ppi)$name), lookup)[, gene]

# Communities ------------------------------
# Community detection based on edge betweenness (Newman-Girvan)
# pdf('output/pdf/graph.pdf', height=10, width=20)

par(mfrow=c(1, 1), mar=c(1,3,1,1), cex=0.6)

set.seed(1)
ceb = cluster_edge_betweenness(hit_ppi) 
V(hit_ppi)$community = ceb$membership
cols = ins(cbPalette, 'white', c(5, 8:9, 11:14))

#svg('output/fig/SFig_ppi_cluster.svg', width=20, height=20)
#dendPlot(ceb, mode="hclust", palette = myCols, cex=1.5, cex.axis=1.5)
plot(hit_ppi, vertex.color=cols[V(hit_ppi)$community],
  vertex.label.cex=2, edge.width=4, vertex.label.family="Helvetica",
  vertex.label.dist=0,  vertex.size=16,
  vertex.frame.color='black', vertex.label.color='black'
)
#dev.off()

The major graph communities, using (Newman-Girvan’s edge betweenness):

  1. The light blue nodes are central and seem to be enriched in regulatory genes.
  2. The grey nodes seem to be a signal transduction from membrane to the regulatory genes (in light blue).
  3. Orange and pink nodes are mostly related to metabolism; those in orange have more membrane related terms than the pink nodes.

Screen Validation

We validate the sensitivity of identified candidate mutants in low throughput and in liquid LB medium at pH 7.4. Instead of factorial (i.e. checkerboard), we use a fixed ratio design Tallarida et al 1997.

dat = fread('input/dat/raw/pa14_mut_pmb_rif.csv') %>% 
  addDrugRatio() %>% 
  addAUC() %>% 
  addFitness() %>% 
  remDissimilarExperimentalConditions() %>% 
  addPA14MutAnnotation()

dat[, dose:= ifelse(cond_f %in% c('PMB', 'combo'), d1, d2)]
dat_lst  = split(dat, dat$mut)
fit = lapply(dat_lst, fitPA14MutDR, curveid = cond_f)

# Loewe response surface ------------------------------
# rs = getPa14RS(dat)

# We'll use a snapshot of results from 2021-01-02
rsl = readRDS('input/dat/rds/2020-01-02_rs_pa14.rds') %>% lapply(., function(x) x$rsl)
foo = lapply(rsl, function(x) summary(x$maxR)$totals) %>% rbindlist(idcol='mut')
# Mutant numbering is potential source of confusion; hence I assign here
# manually: wild-type is coded as mut 0, fitting of mut 28 failed
foo[, mut := c(0:27, 29:45)]  
#arrange(foo, Syn) %>% print(20)
pltNSynergiesOfValidatedMutants()
# ggsave('output/fig/validated_mutants_synergy_count.svg', w=5, h=10)

Dose-response

#svg('output/fig/SFig_45PA14MutDoseResponses.svg', 12, 16)
#plt45PA14MutantDR()  # dose-responses
#dev.off()

Compare to Loewe’s null

#svg('output/fig/SFig_45PA14MutCompLoeweNull.svg', 12, 16)
# We'll use a snapshot of results from 2021-01-02
#plt45PA14MutantComparisonToLoeweNull('2020-01-02_rs_pa14')
#dev.off()

Table with E. coli orthologs

getValidHitTable() %>% addEcOrthologs() %>% .[,-'locus_pa'] %>% 
  knitr::kable()
PA gene EC ortholog Synergy Description
02150 none Serine phosphatase RsbU, regulator of sigma subunit
relA relA none (p)ppGpp synthetase
26590 none GntR family transcriptional regulator
66480 less Predicted ATPase involved in chromosome partitioning
56840 less Predicted thiol oxidoreductase
lldP glcA less L-lactate permease
60490 less cytochrome c
03760 more sodium:solute symporter
43270 selU more tRNA 2-selenouridine synthase
51310 more Predicted redox protein, regulator of disulfide bond formation
maiA more maleylacetoacetate isomerase
43670 more sensor/response regulator hybrid
pmbA pmbA more PmbA protein
62230 more Predicted kinase
  # knitr::kable(format='latex', booktabs=T) %>% 
  # kableExtra::save_kable('output/fig/table_valid_hit.png')

GO term table

getGoTermTable() %>% 
  knitr::kable() 
PA gene Location Process Function
02150 membrane signal transduction catalytic
03760 membrane transmembrane transport transmembrane transporter
26590 NA regulation of transcription; biosynthesis catalytic; transcription factor
43270 NA tRNA seleno-modification transferase for selenium-containing groups
43670 membrane signal transduction; phosphorylation phosphorelay sensor kinase
56840 NA NA electron transfer; heme binding
60490 NA NA electron transfer; heme binding
lldP membrane lactate transport lactate transmembrane transport
maiA cytoplasm aromatic amino acid metabolism catalytic; protein binding
pmbA NA peptidoglycan biosynthesis; proteolysis metallopeptidase
relA NA ppGpp metabolic process NA
  #knitr::kable(format='latex', booktabs=T) %>% 
  #kableExtra::save_kable('output/fig/table_valid_hit_go.png', font='Verdana')

Selected five mutants

3D

#rs = fread('input/dat/raw/pa14_mut_pmb_rif_cfu.csv') %>% 
#  prepForRS()  %>% 
#  setDT()  %>% 
#  addMutGeneNames() %>% 
#  getRS(., 'rs_6mut')
rs = readRDS('input/dat/tmp/rs_6mut.rds')
rsl = lapply(rs, function(x) lapply(x, function(y) y$rsl))

# lapply(names(rs), sav3D, xcap=10)

Membrane permeability

memb = readRDS('input/dat/rds/permeability.rds')
pltMemb(group_by(memb, PmbConcMgL, Strain))
#ggsave('./output/fig/PA14MembranePermeability.svg', w=15, h=5)
memb = readRDS('input/dat/rds/permeability.rds')
om = addAovDunnetMembrane(memb) 
im = addAovDunnetMembrane(memb, membrane = 'IM') 
getDunnetResultsMemb(om, lev = 1)
## $PmbConcMgL
## [1] 0.1
## 
## $RifConcMgL
## [1] 0
## 
## $Dunnett
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`wild-type`
##                      diff     lwr.ci    upr.ci   pval    
## 66480-wild-type -2.106667 -6.2881962 2.0748629 0.4809    
## 02150-wild-type -3.770000 -7.9515296 0.4115296 0.0819 .  
## relA-wild-type  -1.123333 -5.3048629 3.0581962 0.8942    
## 26590-wild-type -1.423333 -5.6048629 2.7581962 0.7821    
## 43270-wild-type  3.376667 -0.8048629 7.5581962 0.1304    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## $RobustDunnett
##  contrast            estimate   SE df t.ratio p.value
##  66480 - (wild-type)    -2.11 1.92 12  -1.098  0.6889
##  02150 - (wild-type)    -3.77 1.91 12  -1.974  0.2417
##  relA - (wild-type)     -1.12 3.01 12  -0.373  0.9773
##  26590 - (wild-type)    -1.42 1.92 12  -0.743  0.8698
##  43270 - (wild-type)     3.38 1.92 12   1.758  0.3293
## 
## P value adjustment: dunnettx method for 5 tests
getDunnetResultsMemb(im, lev = 1)
## $PmbConcMgL
## [1] 2
## 
## $RifConcMgL
## [1] 0
## 
## $Dunnett
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`wild-type`
##                       diff     lwr.ci      upr.ci    pval    
## 66480-wild-type  -7.536667 -16.605398   1.5320642  0.1162    
## 02150-wild-type  -9.753333 -18.822064  -0.6846025  0.0341 *  
## relA-wild-type  -25.990000 -35.058731 -16.9212692 2.2e-05 ***
## 26590-wild-type -29.013333 -38.082064 -19.9446025 5.3e-07 ***
## 43270-wild-type  10.973333   1.904602  20.0420642  0.0169 *  
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## $RobustDunnett
##  contrast            estimate   SE df t.ratio p.value
##  66480 - (wild-type)    -7.54 2.79 12  -2.704  0.0735
##  02150 - (wild-type)    -9.75 5.71 12  -1.708  0.3523
##  relA - (wild-type)    -25.99 3.07 12  -8.471  <.0001
##  26590 - (wild-type)   -29.01 2.14 12 -13.588  <.0001
##  43270 - (wild-type)    10.97 2.72 12   4.041  0.0069
## 
## P value adjustment: dunnettx method for 5 tests

Supplementary Membrane Assays

PMB resistant isolates

dat = 'input/dat/raw/pmb_resistant_isolates.csv' %>%  
  fread() %>% 
  setnames(., 'strain', 'mut') %>% 
  prepForRS(swapnames=T) %>% 
  setnames(., 'mut', 'strain') %>%
  setDT()

# rs = getRS(dat, 'rs_pmbR_intra')
# lapply(names(rs), sav3D)

Change in EC90 of RIF

dat = loadDdiData('input/dat/rds/ddi_all_strains.rds')

dat_fit = anti_join(dat, dat[c1==0 & c2!=0])[, .(data = list(.SD)), .(strain, cond)]

dat_fit[, fit := 
  map(data, function(x) {try(
    drm(effect ~ d1, d2, data = x, 
    fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
    lowerl = c(-4.5, NA, NA))
  )})
]
## Error in optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits,  : 
##   non-finite value supplied by optim
## Error in drmOpt(opfct, opdfct1, startVecSc, optMethod, constrained, warnVal,  : 
##   Convergence failed
idx = sapply(dat_fit$fit, function(x) class(x) != 'try-error')

dat_fit[idx, ed90  := map(fit, function(x) list(try(ED(x, 90, interval='fls', display=F) %>% 
      round(3))))]

dat_fit[idx, par90 := lapply(unlist(dat_fit[['ed90']], recursive=F), function(x) try(as.data.table(x, keep.rownames='level')))]

par90 = lapply(1:nrow(dat_fit[idx, ]), 
  function(x) dat_fit[idx, ][['par90']][[x]]) %>% 
  rbindlist(idcol='id') %>% 
  merge(dat_fit[, .(strain, cond)][idx][, id := 1:26], .) %>% 
  lowerCaseNamesDT() %>% 
  .[, level := gsub('e:|:90', '', level)] %>% 
  .[, effect:= 90]

set.seed(1)
res = par90 %>% 
  group_by(id) %>% 
  mutate(rel_IC = estimate/head(estimate, 1),
    rel_IC = ifelse(rel_IC > 2, 2, rel_IC)
    ) %>% 
  mutate(level = ifelse(level < 1 & level != 0, 0.1, level)) %>% 
  filter(as.numeric(level) <= 1) 

# res %>% group_by(level) %>% summarise(ec90=median(estimate), quantile(estimate, 0.2), ci=gci(estimate), lwr=ec90/ci, upr=ec90*ci) %>% dplyr::select(-ci)
#   level  ec90   lwr   upr
# 1 0      7.41 6.34   8.65
# 2 0.1    3.20 2.32   4.41
# 3 1      2.18 0.966  4.92
# res %>% filter(level==0.1) %>% pull(estimate) %>% quantile(0.2) 
#   20% 
# 0.939 

ggplot(res, aes(level, estimate)) + 
  geom_jitter(width=0.2, pch=21, aes(fill=cond), col='grey70') +
  geom_boxplot(alpha=0.5, outlier.shape=NA) +
  scale_y_continuous(name=expression('EC'[90]^'RIF'*', mg/L'), 
    trans='log2', breaks=c(16, 8, 4, 2, 1, 0.5, 0.25)) +
  scale_x_discrete(name='PMB, xMIC', labels=c('0', '< 1', '1')) +
  coord_cartesian(ylim = c(0.2, 16)) +
  labs(fill = '') +
  theme(aspect.ratio=1)
#ggsave('output/fig/rif_ec90_upon_pmb.svg', width=5, height=4, bg='white')